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Abstract Cranial endocasts can be used as a reliable 
proxy for brain size, reflecting the evolutionary and 
environmental selection pressure of species. Although 
studies on endocasts in amphibians have increased in 
recent years, those performed on endocasts of Anura 
are comparatively rare, especially at the intraspecific 
level. Here, using a high-altitude endemic toad— 
Scutiger boulengeri—as a model, through the application 
of integrative methods (morphology, anatomy, 
phylogeny, and ecology), we studied intraspecific 
variations in endocast morphology and explored its 
driving forces. Three-dimensional reconstruction and 
the brain-to-endocranial cavity (BEC) index suggested 
that the endocast of S. boulengeri can reflect brain 
morphology to a large extent. Elliptic Fourier analysis 
and principal component analysis revealed great 
variability in the cranial endocast morphology among 
individuals, as well as the variation concentrated in 
the regions of telencephalon and optic tectum. In the 
species, individuals with large bodies are accompanied 
by a larger endocast size; the relative endocast sizes 
have significant clade differences but no sexual 
dimorphism. Additionally, the relative endocast sizes 
of S. boulengeri were not associated with phylogenetic 
history and aquatic preference but were positively 
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correlated with altitude and negatively correlated with 
oxygen content, temperature, and precipitation factors 
(annual mean temperature, temperature seasonality, 
annual precipitation, and precipitation seasonality). 
These findings suggested that high-altitude and 
extreme environmental conditions acted as important 
selective forces in morphological variation of the 
cranial endocast of S. boulengeri. 


Keywords aquatic preference, cranial endocast, ecological 
adaptation, phylogenetic history, shape variation 


1. Introduction 


Cranial endocasts, usually defined by the cranial skeleton 
and a container of the brain, share meningeal tissue with the 
brain, probably reflecting varying levels of neuroanatomical 
resolution for different brain regions (Balanoff and Bever, 
2020). The study of cranial endocast morphology contributes 
to conveying important morphological details (shape, extent, 
size, and distribution) for cranial tissues and advances the 
understanding of the evolutionary potential of the brain 
as well as its constraints (Balanoff et al. 2010; Allemand et 
al., 2017). Endocasts serve as empirical windows into the 
history of a phylogenetic stem and play an important role in 
constructing and testing evolutionary hypotheses, informing 
the evolutionary history of neuroanatomical structure and its 
potentially complex evolutionary origin (Marcucio et al., 2011). 
Studies have concluded that the cranial endocast (endocast) does 
provide a reliable morphology proxy for the brain in birds and 


Asian 
Herpetological 
Research 


270 


mammals (brains largely fill the endocranial space with a tight 
environmental correspondence with their braincases) (Iwaniuk 
and Nelson, 2002; Finarelli, 2006). However, in poikilotherms, 
the relationship between the brain and endocast size is much 
more complex (Olori, 2010; Clement et al., 2021). Regardless of 
fish, amphibians, or reptiles, they all present a highly variable 
BEC (brain-to-endocranial cavity) index between groups or 
taxa (Jerison, 2009; Balanoff and Bever, 2020). In amphibians, 
frogs generally have relatively larger brains than salamanders 
(Northcutt, 2002). Recent studies have shown that caecilians 
(68%-74%) and anurans (49%-63%) had brains that filled their 
endocasts to a greater extent than salamanders (38%-41%) 
(Clement et al., 2021), which means that endocasts can be a 
better response to brain morphology in anurans. Clearly, 
however, this conclusion is based on limited data, and brain- 
endocast relationships still need to be further assessed in 
different groups or taxa of amphibians. 

As the central system of learning, memory, cognition, and 
information processing, the brain is an indicator reflecting 
the state of species evolution and development and plays an 
important role in the process of animal adaptation (Sol et 
al, 2005). Variations in vertebrate brain size and structure 
are believed to be related to evolutionary history, ecological 
environment, sexual selection, behavior, developmental 
process, and life history factors (Pollen et al., 2007; Isler and 
Schaik, 2012; Fitzpatrick et al, 2012; Mai et al., 2020; Yu et al, 
2018; Liao et al., 2015). Studies on the adaptive evolution of 
brain variation mainly focus on the interspecific level and 
rarely involve intraspecific level (Gonda et al., 2013; Axelrod 
et al. 2018). Similarly, phylogenetic and ecological signals can 
explain the variations of endocast compared to brain evolution, 
however, only a few studies have focused on the relationship 
between endocast morphology and phylogenetic history or 
environments in macroevolution of ectotherm (Allemand et al., 
2017; Allemand et al., 2019). 

The Qinghai-Tibet Plateau (QTP) is “the third pole of the 
world”, whose average elevation exceeds 4000 m, and QTP 
is one of the important global biodiversity hotspots (Myers 
et al, 2000; Favre et al, 2014). The uplift of the Qinghai-Tibet 
Plateau has caused dramatic changes in climate and ecology, 
which has had a significant impact on the plateau ecosystem 
and has provided a unique opportunity to study the process 
of speciation and ecological adaptation (Favre et al., 2014; 
Wen et al, 2014). To adapt to extreme environments at high 
altitudes, amphibians living in the QTP have evolved a series 
of adaptive phenotypes and physiology, which means that 
these amphibians have become an ideal animal group to study 
morphological adaptation and evolution. 

The Xizang Alpine Toad Scutiger boulengeri (Bedriaga, 1898) 
is a member of the family Megophryidae with an elevation 
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range of 2100-5270 m on the Qinghai-Tibet Plateau (Song et 
al, 1990; Fei et al., 2009; Subba et al., 2015). As one of the most 
widely distributed species with the largest elevation span 
on the Qinghai-Tibet Plateau, information on the species 
morphological, ecological, and phylogenetic information of S. 
boulengeri has received increasing attention (Fei and Ye, 1987; 
Ye et al., 1992; Li et al., 2009; Chen et al., 2009; Hofmann et al., 
2017; Lin et al, 2021). Fei et al. (1987) qualitatively analyzed the 
differences in the webbing and skeleton of the Xizang Alpine 
Toad from different distributions and suggested that during the 
process of adapting to a high elevation and low temperature 
climate, some species gradually adapt to burrowing or aquatic 
environments. Li et al. (2009) indicated that rivers, as an effective 
barrier, led to the differentiation of the internal clades of the 
species. Lin et al. (2021) demonstrated that niche divergence in 
S. boulengeri is caused by niche expansion accompanied by key 
morphological innovations. An extreme environment drove 
sensory adaptations in auditory, visual, and other sensory 
systems, which may be accompanied by changes in the central 
nervous system, including shifts in overall brain size as well as 
modifications in the size and shape of structures of the brain 
(Beltz, 2019). Considering that the endocast can be used as a 
proxy for the brain, herein, to better understand amphibian 
adaptation and evolution to extreme environments at high 
altitudes, using morphological (3D reconstruction and geometric 
morphometry), anatomical, phylogenetic, and ecological 
methods, we first present a quantitative and integrative 
analysis of the endocast morphological variation of S. boulengeri 
and insight into the driving forces from multiple-dimensional 
aspects (evolutionary history, ecological environment, sexual 
selection, etc.). 

Our main goals are to investigate whether the shape and size 
of endocast in S. boulengeri reflect brain morphology and explore 
the morphological variation of the endocast and its possible 
influencing factors. First, we checked the relationship between 
cranial endocasts and brain morphology and quantified 
shape and size variation. Then, we compared the endocast 
shape and size variation among different genders and clades. 
Finally, to explore its driving forces, we tested the correlations 
of phylogenetic history, sexual selection, and environmental 
factors (altitude, latitude, oxygen, habitat preference, and 
climate variable) with endocast size variation. 


2. Materials and Methods 


2.1. Sample preparation All specimens in this study are from 
the Amphibians and Reptiles in the Herpetological Museum 
of the Chengdu Institute of Biology, Chinese Academy of 
Sciences (CIB, CAS) and all procedures followed the guidelines 
established by the Animal Care and Use Committee of the 
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Chengdu Institute of Biology, Chinese Academy of Sciences 
(Permit No.: 2017-AR-JJP-03). The 13 collection sites cover 
three provinces and one autonomous region in the Qinghai- 
Tibet Plateau and its surrounding areas in China. We collected 
altitude, longitude, and latitude for each site (Table 1). The sex 
of the specimens was identified by the spiny glands in the chest 
and fingers (males), unlaid eggs in the abdomen (females), and 
limb shape according to Fei et al. (2009). All specimens were 
fixed in formalin, and we finally selected 98 individuals (67 
males and 31 females) for this study. 


2.2. Measurements and 3D reconstruction The body size 
(SVL: snout-vent length) of specimens (67 males and 31 females) 
was measured by a digital Vernier caliper, accurate to 0.01 mm. 
We scanned 64 individuals (48 males and 16 females) for 
endocast variation analysis. A micro-CT machine (PerkinElmer, 
Quantum GX, in CIB Herpetology Laboratory) was used to 
scan the skull at 70 kV, 88 pA, 14 min per individual. Three- 
dimensional (3D) reconstruction and image processing of the 
endocast were performed using MIMICS v. 21 (Materialise 
Medical Co., Belgium). Each lateral view and dorsal view of the 
endocast was exported from Mimics v. 21 at the same degree 
in “bmp” format, respectively. CL (maximum length: distance 
between the middle anterior-most part of the olfactory bulbs 
and the middle posterior end of the endocast), CW (maximum 
width of the endocast), and CH (maximal height of the 
endocast) of the endocast were measured with the “Measurement 
tool’, and the scanning volume of the endocast (3D volume of 
endocast, V3d) was generated using Mimics v. 21. The mean of 
three measurements was taken, accurate to 0.01 mm. 

To check the brain-endocast relationship, nine CT scanned 
samples (7 males and 2 females) from the site of Tibet 
Autonomous Region and Sichuan Province were selected for 
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skull dissection to obtain the entire brain tissue. The dissected 
brain volume (Vim) was measured using the immersion 
method (Scherle, 1970). The immersion method operation steps: 
The detached brain was placed in an empty tube, and a fixed 
volume of ultrapure water was added with a pipettor until the 
tube was full (total capacity 200 pL). The volume of the brain 
was obtained by subtracting the volume of the water added 
from the volume of the tube (accurate to 1 pL). Measurements 
were repeated three times and averaged to increase the precision 
of the measurement. To negate interobserver variability, all 
measurements were conducted by the same person. 


2.3. Phylogenetic reconstruction We reconstructed the 
intraspecific phylogenetic relationship of S. boulengeri based 
on 16S and COI mitochondrial gene fragments (sampling sites 
are shown in Table 1). The GenBank accession numbers are 
OK544538-48 (16S), OL589209-10 (16S), OK584750-60 (COI) 
and OL597983-84 (COI). We aligned the sequences using 
multisequence alignment in MEGA v.7 (Kumar et al., 2016). 
Before constructing a phylogenetic tree, PartitionFinder v.2.1.1 
(Lanfear et al., 2016) was used to select the best nucleotide 
substitution model for the 16S and CO1 sequences. The optimal 
model was selected based on the Bayesian information criterion 
(BIC) and the greedy algorithm. The maximum likelihood 
method (ML) and Bayesian inference (BI) method were used to 
construct gene trees to analyze the intraspecific phylogenetic 
relationship of species. The BI tree was constructed using 
MrBayes v.3.2.2 (Ronquist et al, 2012) with four Markov chain 
Monte Carlo (MCMC) runs of 10 million generations, sampling 
every 1000 generations and discarding the first 25% as burning. 
ML analysis was performed in RAxML v.7.0.4 software 
(Silvestro and Michalak 2011) with 1000 bootstrap replicates for 
branch confidence evaluation. 


Table 1 Information for sampling locations and sample size of S. boulengeri. Nine complete anatomical brain samples were selected from bold 


sites. 


Regions Sites Sample size Longitude (°E) Latitude (N) Altitude (m) Sites abbreviations 
Sichuan Pengbuxi town, Kangding county 9 101.52860 29.80232 3374 KDPBX 
Sichuan Shade town, Kangding county 5 101.53107 29.49257 3789 KDLB 
Sichuan Sangdui town, Daocheng county 5 100.08412 29.26589 4052 DC 
Tibet Rongbusi, Dingri county 5 87.33460 28.82875 4849 DR 
Gansu Xinglongshan, Yuzhong county 7 104.01013 35.48745 2383 GSYZ 
Sichuan Qialangduoji, Muli county 6 100.51840 28.35921 4359 ML 
Qinghai Tanggu town, Tongde county 5 100.36581 35.35801 3404 QHTD 
Qinghai Shanglaxiu town, Yushu country 5 96.64868 32.77412 4055 QHYS 
Tibet Lawushan, Mangkang county 4 98.54653 29.94484 4187 MK1 
Tibet Luoni town, Mangkang county 3 98.51143 29.69661 4350 MK2 
Tibet Jida town, Basu county 5 96.70958 29.74490 4186 BS 
Sichuan Zhusang town, Yajiang county 1 101.37644 30.03201 4292 YJ1 
Sichuan Decha town, Yajiang county 4 100.74734 29.70842 4028 YJ2 
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The estimation of divergence times of the species was carried 
out in BEAST version 1.7.5 (Suchard et al., 2018). The “uncorrected 
lognormal relaxed clock” model and Yule speciation model 
for branching rates were set in BEAUti. Early diversification 
between Oreolalax and Scutiger occurred during the early 
Eocene (the most recent common ancestor of Oreolalax and 
Scutiger was approximately 53 Ma) (Hofmann et al, 2017), which 
was used as a calibration point to estimate the intraspecific 
clade divergence in S. boulengeri. We performed a Markov chain 
Monte Carlo simulation run for 10 million generations, and 
trees were sampled every 1000 generations. For each of the tree 
statistics in the program Tracer version 1.7.1 (Rambaut et al., 
2018), the effective sample size (ESS > 200) indicated satisfying 
convergence of the Bayesian chain and adequate model mixing, 
A maximum clade credibility (MCC) tree was generated with 
mean node heights and a 10% burn-in through TreeAnnotator 
version 1.7.5.(Drummond et al., 2012). 


2.4. Environment data processing Data on environmental 
factors included habitat factors (quantified by the aquatic 
preference index), altitude, latitude, oxygen content, UV-B 
intensity, and bioclimatic factors (Table S1). Studies have 
suggested that web degree can reflect aquatic preferences in 
amphibians and reptiles (Wei et al., 2009; Foth et al., 2019). By 
default, the larger the web ratio/degree is, the stronger the 
aquatic preference. We quantified aquatic habitat preferences 
using ImageJ (Abramoff et al., 2004) to measure the area and 
calculate the ratio of the webbed area between the third and 
fourth toes. The measurement method of the proportion of 
webbed area was slightly modified based on Yang (2018) (Figure 
S1). The altitude information was collected from the latitude 
and longitude of species sites and the extracted oxygen content 
in the atmosphere (1980-2019) from the National Tibetan 
Plateau Data Center (Xin et al, 2021). Three temperature factors, 
annual mean temperature (bio01), mean monthly temperature 
range (bio02), and temperature seasonality (standard deviation 
*100) (bio04); two precipitation factors, annual precipitation 
(biol2) and precipitation seasonality (coefficient of variation) 
(biol15); and two UV-B factors, annual mean UV-B (UVBI1) 
and UV-B seasonality (UVB2), were selected as representative 
climate factors (Table S2). The climate variables included 19 
bioclimate factors from the WorldClim database for current 
climate (1970-2000) (Hijmans et al., 2005) and six UV-B (UVBI- 
UVBB6) variables from the glUV database (Beckmann et al., 
2014). The bioclimate variables have a resolution of 30 arc- 
seconds (~ 1 km at the equator) for each environment layer 
(Fick and Hijmans, 2017), and the UV-B variables have a spatial 
resolution of 15 arc-minutes. To reduce collinearity and overlap 
between climate variables, Spearman correlation analysis was 
used to remove the highly correlated variables (correlation 
coefficient: |r| > 0.7, P < 0.05). 
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2.5. Shape variation analysis Elliptical Fourier analysis 
(EFA) is a suitable method to analyze the irregular shape of 
objects and can provide an accurate morphological feature 
of the individual through the quantification of the profile, 
facilitating polymorphism analysis in the plants and animals 
(Haines and Crampton, 2000; Raveloson et al, 2005). Given that 
the cranial endocast is irregularly shaped and lacks explicitly 
defined homologous landmarks, EFA was used to quantify 
morphological variation in this structure. Shape analyses were 
performed with normalized elliptic Fourier analyses. First, 
cranial endocast was converted to binary, polished objects and 
chain code using the program SHAPE v.13 (Iwata and Ukai, 
2002). Then, CHC2NEF is used to create Fourier harmonic 
values, and we used 30 harmonics to reconstruct the shape. 
This number of harmonics was enough to capture the shape 
of the outline. In SHAPE, principal component analysis of the 
EFD coefficients was performed using the variance-covariance 
matrix of the coefficients, and the software returned the 
PCA scores for each object. The shape variations that can be 
accounted by each PCA were obtained and graphed using the 
“PrinPrint” tool from SHAPE. 

The gender and clade differences of endocast shape were 
performed by stepwise discriminant analysis (SDA) using SPSS 
v. 24 software (SPSS, USA). SDA was performed on the EFDS 
coefficients after endocast shape standardization (coefficients 
without difference were unsuitable for entering into the 
equation and would be excluded) to find linear combinations of 
the original variables that could give best separation between 
groups. The significance of the discriminant function was 
evaluated from Wilks’ lambda statistic (level of significance set 
at 0.05) (Klecka, 1980). 


2.6. Phylogenetic signal testing The species is widely 
distributed on the Qinghai-Tibet Plateau, with some 
divergences in morphology and phylogeny (Li et al., 2009; 
Hofmann et al., 2017; Lin et al., 2021), so we used phylogenetic 
signals to investigate the similarity of traits among clades 
(Blomberg et al., 2003). We estimated the intensity of the 
phylogenetic signal in SVL, web degree, and absolute and 
relative endocast size (residuals from the endocast to SVL) 
using the package “phytools” in R software (R Core Team, 
2014). Phylogenetic signals are the tendency for closely related 
species to display similar trait values as a consequence of their 
phylogenetic proximity (Blomberg et al., 2003). Phylogenetic 
signals calculated by Blomberg’s K-statistic method can measure 
the phylogenetic signal intensity of functional traits and detect 
the correlation between functional traits and phylogenetic 
relationships among species (K values closer to zero correspond 
to a random or convergent pattern of evolution, K < 1 implies 
that there is less phylogenetic signal than expected under 
Brownian motion, while K > 1 indicates a strong phylogenetic 
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signal and conservatism of traits) (Blomberg et al., 2003; Adams, 
2014). 


2.7. Difference analysis and correlation tests First, we test 
the assumption of normality and homogeneity of variances of 
the data. To examine the gender and clade differences in SVL, 
web degree, and cranial endocast, we used two-way analysis 
of variance (two-way ANOVA) and analysis of covariance 
(ANCOVA) and a nonparametric factorial Kruskal-Wallis 
(KW) sum-rank test. Then, Fisher's least significant difference 
and Dunn’s test analysis were used for multiple comparisons 
of the above data. The dimorphism index is directly calculated 
as the ratio of male size to female size minus 1, with a positive 
value indicating larger males and a negative value indicating 
larger females (SDI = (male size/female size) — 1) (Lovich and 
Gibbons, 1992). A general linear regression model (GLM) 
was used to investigate the relationship between immersion 
brain volume (Vim) and V3d, and a regression equation 
was established to verify the effect of scanning volume. To 
examine the correlation between the web degree, altitude, 
latitude, oxygen content, climate factors and endocast sizes, 
we constructed linear mixed models (LMMs) in the R 
package “Ime4’. To compare the relative size, GLM was used 
to investigate the relationship between endocast sizes and 
SVL, and SVL was set as a covariable to control its impact in 
exploring the relationship between endocast sizes and all factors. 
The data (SVL, web degree, and cranial endocast measurements) 
meet the assumption of normality and homogeneity of 
variances, except for the CW data (homogeneity of variance is 
not satisfied in different clades). For group differences in SVL 
and web degree, gender, clades, and their interaction were fixed 
factors. For endocast sizes, gender, clades, and their interaction 
as the fixed factor and the SVL as a covariable control the body 
size effect. For the gender difference of CW data, gender was the 
fixed factor, and SVL was the covariable. The KW sum-rank 
test was used to analyze the difference in CW data size (residual 
data corrected by SVL) between clades. For all GLM and 
LMMs, we ran model diagnostics to test model assumptions, 
including normality and homoscedasticity. The applied LMMs 
are as follows: 

1. Aquatic preference: an LMM involving web ratio as the 
response variable, gender and altitude as fixed variables, with 
clade as a random factor. 

2. Habitat preference: LMMs involving V3d, CL, CH, CW as the 
response variable; web degree as a fixed factor, with clade as a 
random factor. 

3. Altitude and latitude: LMMs involving V3d, CL, CH, and 
CW as the response variables; altitude and latitude as a fixed 
factor, respectively, and SVL was considered a covariate to 
control its impact, with clade as a random factor. 

4. Oxygen content: LMMs involving V3d, CL, CH, CW as the 
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response variable; oxygen content as a fixed factor; SVL was 
considered a covariate to control its impact, with clade as a 
random factor. 

5. Climate factors: Before the correlation analysis between 
endocast sizes and climate factors, GLMs were used to test 
the correlation between altitude, latitude and climate factors, 
respectively. LMMs involving V3d, CL, CH, CW as the response 
variable; climate factors (bio01, bio02, bio04, biol2, biol5, UVBI, 
UVB2) as fixed variables; SVL was considered a covariate to 
control its impact, with clade as a random factor. 


3. Results 


3.1. Phylogeny After manual correction, alignment, and 
trimming, we finally obtained 1125 bp of two mitochondrial 
gene fragments (16S: 502 bp, COL 623 bp). The best substitution 
model for both gene fragments was GTR +G. Based on the 
BI and ML phylogenetic trees, we divided S. boulengeri into 
five clades, including GS, QH, S1 (including the sites from YJl, 
KDPBX, KDLB), S2 (including the sites from ML, DC, YJ2), 
and XQ (including the sites from DR, BS, MK1, MK2, QHYS). 
The resultant topology was consistent with the intraspecific 
phylogeny of S. boulengeri obtained from the analysis of the 
cytochrome b gene (Lin et al, 2021). Divergence times indicated 
that the clades of S. boulengeri diverged during the Miocene 
approximately 13.66 Ma (Figure 1). 


3.2. Morphological comparison between brain and 
endocast In S. boulengeri, the external surface of the brain 
and inner surface of the endocranial cavity are consistently 
buffered by meningeal tissue (dura, pia, arachnoid). We 
obtained different brain regions, including the olfactory 
nerves, olfactory bulbs, telencephalon, pituitary infundibulum, 
diencephalon, optic tectum, cerebellum, and medulla oblongata 
(Figure 2D). A pair of corresponding olfactory nerves extend 
forward, shape elongated; the bulbs olfactory is cylindrical, 
short, transversally undivided; and the telencephalon consists 
of a pair of evaginated cerebral hemispheres (both hemispheres 
are separated by a longitudinal medial sulcus), drop-shaped 
and globoid; the diencephalon is pentagonal from dorsal view, 
partially covered by the posterior edge of the hemispheres. The 
pituitary infundibulum is relatively flat and located on the 
ventral side of the diencephalon. In dorsal view, the optic tectum 
exhibits a pair of domes separated by a median sulcus, nearly 
round in shape, and the cerebellum is relatively narrow and 
thin. The medulla oblongata is sunken-shaped and posteriorly 
connects with the spinal cord (Figure 2D). Mirco-CT scan and 
3D reconstruction showed that the skull of the S. boulengeri 
specimen could be reconstructed completely; however, from 
the endocast, it is not possible to distinguish different brain 
regions. The endocast (without considering the cranial nerves or 
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Oreolalax chuanbeiensis 
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Figure 1 The phylogenetic tree of S. boulengeri with a relaxed clock and Yule speciation process. The number on the node label is the 


estimated divergence time. Numbers inside the branches indicate Bayesian posterior probabilities/maximum likelihood bootstrap 
values, respectively. 


A C 
we Bulbus nerve region 
Olfactory bulb region 


Telencephalon region 


Optic tecta region 


Bulbus nerve 
Olfactory bulb 


Telencephalon 


Diencephalon 
Optic tecta 


Cerebellum 
Medulla oblongata 


Figure 2 Skull, cranial endocast and dissected brain of S. boulengeri. A: represents dorsal view of skull; B: represents lateral view of skull; C: 
represents the ventral, lateral, and dorsal view of cranial endocast (left to right), the Height, Length, and Width represent the measurement 
parameters of endocast; D: represents the ventral, lateral, and dorsal view of brain (left to right). 
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the inner ear) external surface was relatively smooth, although 
there were some individual differences (Figure 2A, B, C). The 
morphological comparison between anatomical brain and 
reconstructed endocast showed that the regions of bulbs nerve, 
olfactory bulb, telencephalon, pituitary infundibulum and optic 
tectum of endocast show high morphological consistence with 
brain shape (Figure 2C, D). The bulb nerve corresponds to the 
anterior-most structure of the brain endocast, and the largest 
part has a groove (corresponding to a groove in telencephalon) 
in the middle and laterally extended. The external morphology 
of the brain endocast does not allow the olfactory bulb and 
telencephalon to be delimited in dorsal view, and the lateral 
view shows a relatively clear delimitation. In the dorsal 
view, the median sulcus of the optic tectum is not visible. The 
posterior margin from the optic tectum in the endocast covers 
the indistinct cerebellum and medulla oblongata. The pituitary 
infundibulum is marked by a small bulge in ventral and lateral 
views (Figure 2C). 

The immersion method and 3D reconstruction were used 
to measure the size of the brain and endocranial cavity, 
respectively. We obtained and compared the brain size with 
the endocast volumes by CT reconstruction in S. boulengeri and 
found that the BEC index was more than 1/2 (the proportion 
of brain to endocast volume in the measured specimen was 
56.38%-68.14%, the average proportion was 61.29%). General 
linear regression analysis showed that there was a significant 
positive correlation between brain size and endocast volume 
(Vim = 0.504*V3d-+6.4737; P = 0.005, adjusted R’= 0.66). 


3.3. PCA of endocast shape variation Principal component 
analysis (PCA) summarized the cranial endocast shape 
variations of S. boulengeri. In the dorsal view of endocast shape, 
the first five principal components (79.4%) mainly represented 
the variation of the whole endocast (Table 2). The first three 
PC axes explained 36.3%, 20.6%, and 9.9% of the observed shape 
variation. PC1 and PC2 described the variation of telencephalon 
region; with PC1 represented the broaden and shorten 
telencephalon region from the negative to positive direction; 
PC2 showed the narrowing telencephalon region from the 
negative to positive direction. PC3 captured the widened optic 
tectum region variation from the negative to the positive 
direction. PC4 was related to extending the forward olfactory 
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bulb region and the general shape variation of the endocast. 
From the PCA, there were obvious overlaps in gender and 
clades. Figure 3 showed that there is a large variation among 
individuals in the endocast of S. boulengeri. 

In the lateral view of the cranial endocast, the first five 
principal components (74.4%) mainly represented the variation 
of the cranial endocast (Table 2). The first three principal 
components together accounted for approximately 60.0% of 
the lateral shape variation. The first, second, and third PC 
axes explained 33.9%, 15.7% and 10.4% of the lateral observed 
variation, respectively. For PCI, telencephalon and optic tectum 
regions were higher in the positive direction, while the shape 
was flatter in the negative direction. For PC2 and PC3, the 
main shape variation focused on the olfactory bulb and optic 
tectum regions. The negative to positive direction in PC2 was 
characterized by the optic tectum flatter to round and round 
to flat in the olfactory bulb. PC3 captured the variation in 
flatting the olfactory bulb and optic tectum regions from 
negative to positive. PC4 mainly represented the flatting 
olfactory bulb from negative to positive. From the results of 
principal component analysis (many obvious overlaps were 
found in gender and clade), there was a large variation among 
individuals in the cranial endocast of S. boulengeri (Figure 3). The 
discriminant analysis was performed on 117 EFDS coefficients 
after endocast outline standardization. The results indicated 
that the cross-validation rates of gender in the dorsal and lateral 
views of the Fourier coefficient were 65.6% and 76.6% using SDA, 
respectively, suggesting that there were some differences in the 
endocast shape of the gender from the dorsal and lateral views, 
but the discrimination accuracy was low. In the discrimination 
of clade outlines, we could not distinguish the clade differences 
using the endocast outlines because the cross-validation rates 
of the dorsal view and the lateral view were 40.6% and 59.4%, 
respectively. All the discriminant functions were valid (Wilks’ 
lambda significance P < 0.05). The discrimination accuracy of 
endocast shape for both genders and clades was low, which also 
reflected the large variation among individuals. 


3.4. Group difference of traits in S. boulengeri 

Sexual dimorphism The SVL and web degree had significant 
sexual dimorphism (SVL: F = 8.638, P = 0.004, web degree: F = 
10.908, P = 0.001). The dimorphism index results showed that 


Table 2 Eigenvalue, contribution rate and cumulative contribution rate of the first five principal components of elliptic Fourier descriptors. 


Dorsal view Eigenvalue Proportion (%) Cumulative (%) Lateral view Eigenvalue Proportion (%) Cumulative (%) 
Rei 1.65E-03 36325 36:325 IKOS 4.98E-04 33.878 33.878 
PC2 9.35E-04 20.586 56.911 PC2 2.32E-04 15.740 49.619 
BGS 4.50E-04 9.910 66.821 RES 1.53E-04 10.399 60.018 
PC4 3.10E-04 6.814 73.635 PC4 1.18E-04 8.035 68.053 
RES 2.61E-04 DST, 72373 RES 9.37E-05 6.369 74.422 
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Figure 3 Visualization of the variations in endocast shape (dorsal and lateral view) in S. boulengeri. A: The blue and black solid lines 
indicate, respectively, the low and high values along the first four PC axes; B-C: represent the shape variation of different clades in dorsal 
view (first three PCs); D-E: represent the shape variation of dorsal view (first three PCs) in both genders; F-G: represent the shape 
variation of different clades in lateral view (first three PCs); H-I: represent the shape variation of lateral view (first three PCs) in both 
genders. Each endocast shape is reconstructed from coefficients calculated by letting the score for the corresponding principal component 
be equal to its mean or its mean plus or minus two times the standard deviation. The endocasts (dorsal and lateral view) represent the 
shape variation from negative to positive along the first three PC axes. The blue arrows point the variation direction from negative to 


positive. 


females had a slightly larger body size than males (SDI = 
-0.016 < 0), whereas males had a larger degree of webbing than 
females (SDI = 0.240 > 0). In S. boulengeri, we found that 
a large body was accompanied by a large endocranial cavity 
(P = 0.017, adjusted R’ = 0.074, Table $3), so body size was used 
to correct the relative endocast size. No sexual dimorphism was 
found in any of the relative endocast sizes (V3d, CL, CH, CW) 
(P > 0.05). 

Clade difference ANOVA results showed that SVL and web 
degree had significant clade differences (SVL: F = 5.715, P < 0.001; 
web degree: F = 12.753, P < 0.001). The relative V3d and CH 
were significantly different in clades (V3d: F = 8.257, P < 0.001; 
CH: F = 5.312, P = 0.001), while the relative CW and CL were 
not significantly different among clades (P > 0.05). The results 
suggested that endocast sizes in different clades differentiated to 


‘web degree 


some degree (Figure 4). 


3.5. Correlation between phylogenetic history, 
environmental factors and endocast sizes in S. boulengeri 
Phylogenetic signal To analyze whether features evolved 
with the Brownian motion model, we used the average value of 
each site to detect the phylogenetic signals. Body size (SVL), web 
degree and absolute endocast volume had weak phylogenetic 
signals (K < 0.5, P < 0.05). Controlling the body size, the results 
showed that there was no significant phylogenetic signal in 
relative cranial endocast sizes (residual) (Table 3). 

Aquatic preference To quantify the habitat preference 
factors, an indirect index (web degree) was selected to reflect 
the differences in aquatic preference in S. boulengeri. The 
LMM results suggested that there was a significantly positive 
correlation between web degree and altitude (male: t value = 
3.039, P = 0.004; female: t value = 2.857, P = 0.007) (Figure 5A), 
indicating that the higher the altitude, the larger the webbing, 
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Figure 4 The morphological measurements in clades of S. boulengeri. A: SVL, B: Web degree, C: relative 3D Volume of endocast (V3d~ 
residual), D: relative endocast height (CH~residual). 


Table 3 The phylogenetic signal (K-value) of all traits in S. boulengeri. 


Traits K P (K) Traits K P (K) 
V3d-a 0.180 0.041 V3d-r 0.142 0.107 
CL-a 0.088 0.420 CL-r 0.075 0.539 
CW-a 0.145 0.107 CW-r 0.106 0.266 
CH-a 0.104 0.292 CH-r 0.091 0.392 
SVL 0.244 0.015 Web degree 0.164 0.043 


Note: V3d represents 3D volume of endocast; CL represents maximum endocast length; CW represents maximum endocast width; CH 


«oo» 


represents maximum endocast height; SVL represents body size. The trait with “-a” means the absolute size and with the “-r” means the 


relative size. 


suggesting an increased aquatic preference with altitude. 
Meanwhile, we investigated the correlation between aquatic 
preference (web degree) and endocast sizes in the species. 
The results confirmed that the relative endocast sizes had no 
significant correlation with aquatic preference (P > 0.05) (Table 
S4). 

Altitude and Latitude In examining the correlations between 
altitude, latitude and endocast sizes, the results showed that the 
relative V3d and CH had a positive correlation with altitude 
(V3d: P = 0.032; CH: P = 0.022), and the relative CW and CL had 
no correlation with altitude (P > 0.05) (Table S5). No significant 
correlation was found between relative endocast sizes and 
latitude (P > 0.05) (Table S5). The results showed a trend that 


the relative endocast sizes increased with increasing altitude 
(Figure 5B). 

Oxygen content Oxygen content, an associated factor with 
elevation change, was further tested with endocast sizes. The 
results showed that the relative V3d and CL had a negative 
correlation with oxygen content (V3d: P = 0.010; CL: P = 0.034), 
and the relative CW and CH had no correlation with oxygen 
content (P > 0.05) (Table S6). That was, the relative endocast 
volume of the species increased with the gradual reduction of 
oxygen content under the environment of the plateau (Figure 
6). 

Climate factors GLM results showed that altitude was 
negatively correlated with bio01, bio04, biol5 and positively 
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Figure 5 Correlation between web degree, relative endocast volume and altitude (ALT) in S. boulengeri. A: web degree with altitude (ALT) 
in both genders; B: relative endocast volume (V3d~residual) with altitude (ALT). 
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Figure 6 Correlation between oxygen content (OC) with the 
relative endocast volume (V3d~residual) in S. boulengeri. 


with bio02, UVB1, UVB2, respectively. The latitude was 
negatively correlated with biol2, UVB1, UVB2 and positively 
with bio04, respectively (Table S7). Then we analyzed the 
relationship between multiple climate factors (bio01, bio02, 
bio04, biol2, biol5, UVB1, UVB2) and endocast sizes variation. 
The linear mixed-effects models revealed that the relative 


endocast sizes were related to climate factors. In further 


exploring the relationship between the temperature factor 
and relative endocast sizes, we found that relative CH was 
significantly negatively correlated with bio01 (annual mean 
temperature) (P = 0.039, Table S8), while other measurements of 
endocasts were unrelated to the factor (P > 0.05, Table S8). All 
endocast sizes had no correlation with bio02 (mean monthly 
temperature range) (P > 0.05, Table S8). The relative V3d was 
significantly negatively correlated with bio04 (temperature 
seasonality) (P = 0.018) and relative CW showed a weak 
negative correlation with bio04 (P = 0.05), while the relative 
CL, and CH were unrelated to this factor (P = 0.05) (Table S8). 
For the precipitation factors, the results showed that relative 
V3d and CH were negatively correlated with biol2 (annual 
precipitation) and biol5 (precipitation seasonality), respectively 
(Table S8). Moreover, the UV-B factors (annual mean UV-B 
and UV-B seasonality) had no correlation with endocast sizes (P 
> 0.05, Table S8). The influence of temperature and precipitation 
factors on the endocast was evident, and the smaller seasonality 
variation was favorable to generate a larger cranial endocast 
(Figure 7). 


4. Discussion 


4.1. The endocranial cavity can reflect brain morphology 
in S. boulengeri Using CT scanning combined with 3D 
reconstruction, we successfully reconstructed the endocranial 
space of S. boulengeri and obtained cranial endocast morphology. 
In morphological comparisons between the brain and endocast, 
the regions of the olfactory bulb, telencephalon and optic 
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Figure 7 Correlation between climate factors with the relative endocast volume (V3d~residual) and relative endocast height (CH~ 
residual) in S. boulengeri. A: bio04 (Temperature seasonality) with V3d (residual), B: bio15 (precipitation seasonality) with V3d (residual), C: 


bio12 with CH (residual), D: bio15 with CH (residual). 


tectum of the endocast showed high morphological similarity 
with brain shape. Moreover, there was a significant positive 
correlation between the endocast cavity and brain size. 
Therefore, we interpreted that the endocranial cavity of S. 
boulengeri can estimate brain size to a large extent (average BEC 
index is 61.29%, range of 56.4%-68.1%). This is the first report 
on the Mesobatrachia species, and the result is consistent with 
the results of Challands et al. (2020) and Clement et al. (2021) 
on Neobatrachia species (49%-63%). Our study provides new 


evidence for a relatively tight relationship between the endocast 
and brain morphology and refutes the argument of a low BEC 
index in amphibians (Jerison, 2009). 


4.2. Shape variation of endocasts In this study, we applied 3D 
reconstruction and elliptical Fourier analysis to investigate the 
variation in cranial endocasts in S. boulengeri. From the lateral 
and dorsal views, the most distinct variations of individuals 
located in the region of telencephalon and optic tectum. The 


endocast morphology varied greatly among individuals, and it 
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was not very well classified in genders and clades. Our results 
showed no clear pattern in the endocast shape, but the presence 
of differences was undeniable, more samplings and innovative 
comparative quantitative methods (higher-resolution scan or 
dice CT, 3D geometric morphometrics) can more accurately 
reconstruct and quantity the shape of the structures in the 
future. 


4.3. No sexual selection in relative endocast sizes In brain 
morphology studies of birds and mammals, the intensity of 
sexual selection is usually treated as a key component affecting 
the evolution of brain size (Garamszegi et al, 2005; Fitzpatrick 
et al, 2012). In amphibians, fewer studies have investigated how 
sexual selected pressures influence investment in brain tissue 
(Zeng et al, 2016; Mai et al, 2020). Zeng et al. (2016) indicated that 
sexual selection does not play a significant role in the evolution 
of brain size in Anura, although courtship behavior and the 
mating system affect the size evolution of the olfactory brain. 
Mai et al. (2020) found that relative brain size was positively 
correlated with intensity of male mate competition (operational 
sex ratio (OSR), spawning-site density and male forelimb 
muscle mass). Our study suggested that the web degree and 
SVL of S. boulengeri have significant sexual dimorphism (P < 
0.001), while there was no significant sexual dimorphism for 
relative endocast sizes (V3d, CL, CW, CH) (P > 0.05). Therefore, 
we speculated that sexual selection is not an important factor in 
endocast size evolution within the species. 


44, Evolutionary instability of endocast sizes It is common 
to use phylogenetic signals to predict features of species that 
have not been studied, as well as features of their close relatives 
(Adams 2014). Studies have demonstrated that the evolution 
of morphological differences between clades or species could 
be the consequence of evolutionary constraints (Vidal-Garcia 
et al, 2014). Additionally, morphological characteristics are 
considered to be evolutionarily malleable and more susceptible 
to environmental influences or interactions between ecological 
conditions and historical evolutionary context (Axelrod et al., 
2021). S. boulengeri harbors a large amount of genetic variation, 
and the mtDNA or nuclear genes of the species have a strong 
geographical pattern (Li et al, 2009; Chen et al., 2009; Hofmann 
et al, 2017; Lin et al, 2021). Our study identified S. boulengeri into 
five clades according to the phylogenetic relationship (Figure 
1). We examined the relationship between phylogenetic history 
and traits (SVL, web degree, cranial endocast) of the species. 
The results showed that the SVL, web degree and absolute V3d 
had weak phylogenetic signals, and the relative endocast sizes 
had no phylogenetic signals in S. boulengeri. Lin et al. (2021) noted 
that the SVL in the species has a significant phylogenetic signal, 
and our study is consistent with this result. The weak signal in 
absolute endocast volume and lack of phylogenetic signals in 
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the relative sizes of the endocast suggested that the endocast 
sizes within the species tended to be evolutionarily unstable and 
that the variation could be influenced by other factors, such as 
environmental factors. 


4.5. No effect of aquatic preference on endocast sizes 
Differences in the habitat environment have been proven 
to influence the morphology of the brain, endocast, or skull 
(Taylor et al, 1995; Liao et al, 2015; Allemand et al., 2017; Paluh 
et al, 2020). For instance, arboreal families (or subfamilies) 
had a significantly larger brain than terrestrial or aquatic 
families (or subfamilies) in Anura (Taylor et al., 1995). Paluh et 
al. (2020) suggested that skull shapes differ among frogs that 
occupy distinct microhabitats. Previous studies supposed that 
during the process of adapting to the QTP, Scutiger species 
gradually adapt to burrowing or aquatic environments (Ye et 
al., 1992). Our study found that the web degree of S. boulengeri 
increases significantly with elevation, which implied a more 
aquatic adaptation trend increasing with altitude. Furthermore, 
we investigated the relationship between habitat (aquatic 
preference) and endocast sizes in S. boulengeri; however, the 
statistics indicated that there was no correlation between 
endocast sizes and aquatic preference. Powell and Leal (2014) 
suggested that differences in relevant habitat complexity within 
a given habitat type might not be sufficient to favor divergence 
in the volume of the overall brain or its constituent structures, 
especially in closely related species. Thus, we concluded that 
the complexity of aquatic preference habitats alone might not 
be sufficient to favor variation in endocasts at the intraspecific 
level. 


4.6. Effect of environmental factors on relative endocast 
sizes Generally, in high-altitude areas, lower temperatures, 
less food supplies, shorter breeding seasons, and greater 
climatic unpredictability have a great influence on the 
physiology, behavior, and life history of animals, especially 
in ectotherms (Badyaev and Ghalambor, 2001; Liao and Liu, 
2008). The Qinghai-Tibet Plateau is characterized by extreme 
environments, such as low oxygen, low temperature, less 
precipitation, and strong UV-B radiation, which bring multiple 
pressures and challenges to the local species (Wen, 2014). In 
our study, the integrated influence of multiple environmental 
factors on the variation in endocast size was found. The 
endocast size of S. boulengeri was positively correlated with 
altitude and negatively correlated with oxygen content, annual 
mean temperature, and annual precipitation but not correlated 
with the latitude and UV-B factor. Large brain facilitates 
the construction of behavioral responses to unusual, novel or 
complex socioecological challenges (cognitive buffer hypothesis, 
CBH) (Dunbar, 1998; Sol, 2009; Sol et al, 2007). Amphibians can 
evolve larger brains or brain regions as an adaptation for a 
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better chance of colonizing harsher environments (Amiel et al., 
2011; Huang et al., 2020). Lin et al. (2021) uncovered that niche 
expansion to a new environment was more representative of 
real niche change than niche unfilling in the shifted range of 
S. boulengeri. Therefore, as the CBH predicted, we believed that 
the cognitive advantage of a relatively large endocast (or brain) 
in the species contributes to increased behavioral flexibility 
and food availability to adapt to complex geographical 
environments and extreme climate conditions or to colonize 
new environments and expand the ecological niche in the QTP 
and its adjacent areas. 

Seasonal environment variability is also thought to have 
contributed to the evolution of the brain (van Woerden et 
al., 2010; Luo et al., 2017). The expensive brain framework 
(EBF) predicts that the brain size decreases with seasonality, 
as the most energetically costly organ, relatively smaller 
brains can provide energy benefits during periods of resource 
scarcity (van Woerden et al., 2010). Our study found that 
endocast size was negatively correlated with temperature 
seasonality and precipitation seasonality, but not with UV-B 
seasonality, reflecting that lower seasonal variation in extreme 
environmental conditions is more conducive to the development 
of a larger endocast or brain in S. boulengeri. Therefore, our 
study supports the EBF on account of endocast size of S. 
boulengeri decreases with seasonal changes (temperature and 
precipitation). The climatic seasonality is influenced by latitude 
and altitude (Chown and Klok, 2003; Martin et al., 2009). In the 
Tibetan Plateau and its adjacent area, seasonal variations in 
temperature are more pronounced at lower elevations than at 
higher elevations (Lu et al, 2010). Our results further confirmed 
that altitude and latitude do affect the temperature seasonality 
and precipitation seasonality (see Table S7). The endocast sizes 
of S. boulengeri was positively correlated with altitude but not 
with latitude. Given the distribution (large elevation spans) of 
S. boulengeri, the higher altitude populations/clades within lower 
climatic seasonality environments have a larger endocast size or 
brain. 


5. Conclusions 


Here, we quantified the BEC index and endocast shape 
variation, exploring the relationship between sexual selection, 
phylogenetic history, and environmental pressures on endocast 
size variation in an alpine lazy toad—Scutiger boulengeri. The 
study showed that the BEC index in this Mesobatrachia 
species was more than sixty percent, which is consistent with 
the results from Neobatrachia. EFA and PCA results showed 
great variability in the brain endocast morphology among 
individuals, and the main variation focused on telencephalon 
and optic tectum. In addition, our study indicated that endocast 
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sizes had no relationship to phylogenetic history, sexual 
selection and aquatic preference, but was strongly related to 
environmental conditions (especially altitude, oxygen content, 
temperature, and precipitation). These findings suggested that 
the extreme environmental conditions in the Qinghai-Tibet 
Plateau and its surrounding areas played an important role in 
cranial endocast variation. 
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Appendix 


Figure S1 Schematic view of the foot and web area of S. boulengeri. Green area represents webbing size and the dark blue dotted line 
represents space size between third and fourth toes. Web degree were calculated by web size/space size. 


Table $1 Environment data list (except the habitat preference factor). 


Abbreviation Variable name Citation 


> 


Latitude This study 


3 


Annual mean UV-B Beckmann et al. 2014 


s 


mean UV-B of the highest month Beckmann et al. 2014 


: 


Sum of monthly mean UV-B during the highest quarter Beckmann et al. 2014 


bio01 Annual mean temperature Hijmans et al. 2005 


bio03 Isothermality Hijmans et al. 2005 


bio05 Max temperature of the warmest month Hijmans et al. 2005 


bio07 Temperature annual range Hijmans et al. 2005 


bio09 Mean temperature of the driest quarter Hijmans et al. 2005 


bioll Mean temperature of the coldest quarter Hijmans et al. 2005 


biol3 Precipitation of the wettest month Hijmans et al. 2005 


biol5 Precipitation seasonality Hijmans et al. 2005 


biol7 Precipitation of the driest quarter Hijmans et al. 2005 


biol9 Precipitation of the coldest quarter Hijmans et al. 2005 


Table $2 Spearman correlation tests for climate variables (the variables correlation coefficient: |r| < 0.7 were retained). 


Climate variables UVB2 bio01 bio02 bio04 biol2 biol5 


UVB2 -0.592** 0.623** 0.224 -0.302* 0.121 


bio02 0.052 -0.512** 0.356** 


bio12 0.225 


Note: * represents the significant 0.01 < P < 0.05, ** represents the significant 0.001 < P < 0.01. 


Table S3 General linear regression model results of endocast sizes and body size (SVL). 


Estimate Std. Error t-value P-value Adjusted-R” 
V3d 0.7234 0.295 2.452 0.017 0.074 
CL 0.03899 0.01216 3.206 0.002 0.128 
CW 0.04315 0.0115 3.753 <0.001 0.172 
CH 0.024149 0.007236 3.337 0.001 0.139 


Note: V3d represents 3D volume of endocast; CL represents maximum endocast length; CW represents maximum endocast width; CH 
represents maximum endocast height. The bold text represents the significance P-value (P < 0.05 or P < 0.001). 


Table $4 Linear mixed model results of endocast sizes and web degree. 


Response variable Fixed effect Estimate Std. Error df t-value P-value 
Web degree 2.957 17.784 58.231 0.166 0.869 
pa SVL 0.386 0.304 60.778 1269 0.209 
Web degree 0.497 0.656 6l 0.758 0.451 
Be SVL 0.038 0.012 61 3.061 0.003 
Web degree -0.199 0.439 49.335 -0.452 0.653 
Sy SVL 0.019 0.008 60.252 2.491 0.016 
cw Web degree 0.713 0.666 14.516 IOT 0.301 
SVL 0.042 0.012 29.774 3.411 0.002 


Note: V3d represents 3D volume of endocast; CL represents maximum endocast length; CW represents maximum endocast width; CH 
represents maximum endocast height; SVL represents body size. 


Table S5 Linear mixed model results of endocast sizes and altitude (ALT) and latitude (LAT). 


Response variable Fixed effect Estimate Std. Error df t-value P-value 
ALT 0.00564 0.00235 12.674 2.404 0.032 
[es SVL 0.37724 0.29094 60.993 1.297 0.200 
ALT 0.00011 0.00006 61 1.837 0.071 
fr SVL 0.03569 0.01207 61 2.957 0.004 
ALT 0.00013 0.00004 5.043 3.281 0.022 
ae SVL 0.01956 0.00697 40.5 2.808 0.008 
NEAL 0.00001 0.00006 4.271 0.193 0.856 
Dr SVL 0.04275 0.01186 24.6 3.606 0.001 
LAN -0.9617 0.5649 5129 -1.703 0.148 
ae SVL 0.4178 0.2996 56.751 1.395 0.169 
LAT -0.02547 0.01689 61 -1.508 0.137 
Pe SVL 0.03843 0.01204 61 3.191 0.002 
LAT -0.01712 0.01222 5.964 -1.401 0.211 
mo SVL 0.019804 0.0075 46.771 2.641 0.011 
LAT -0.01209 0.01619 61 -0.746 0.458 
re SVL 0.04289 0.01155 61 3.714 <0.001 


Note: V3d represents 3D volume of endocast; CL represents maximum endocast length; CW represents maximum endocast width; CH 
represents maximum endocast height; SVL represents body size. The bold font represents the significance (P < 0.05 or P < 0.001) of the 
correlation between endocast sizes and ALT. 


Table S6 Linear mixed model results of endocast sizes and oxygen content (OC). 


Response variable Fixed effect Estimate Std. Error df t-value P-value 


SVL 0.466 0.289 60.427 1613 0.112 
SVL 0.022 0.007 32.851 3.104 0.004 
SVL 0.038 0.012 18.35 3.120 0.006 
SVL 0.043 0.012 21.731 3.625 0.002 


Note: V3d represents 3D volume of endocast; CL represents maximum endocast length; CW represents maximum endocast width; CH 
represents maximum endocast height; SVL represents body size. The bold font represents the significance (P < 0.05) of the correlation between 
endocast sizes and OC. 


Table $7 General linear regression model results of the correlation between altitude (ALT), latitude (LAT) and selected climate variables, 
respectively. 


Formula Estimate Std. Error t-value P-value 


ALT ~bio02 171.56 68.54 2.503 0.015 


ALT ~bio12 -0.5248 0.6061 -0.866 0.390 


ALT ~UVB1 0.57937 0.0526 11.020 <0.001 


LAT ~bio01 0.1219 0.1168 1.044 0.300 


LAT ~bio04 0.01651 0.00245 6.740 <0.001 


LAT ~bio15 0.02627 0.02581 1.018 0.313 


LAT ~UVB2 -0.00004 0.00001 -5.043 <0.001 


Note: The bold text represents the significance P-value (P < 0.05 or P < 0.001). 


Table S8 Linear mixed model results of endocast sizes and climate variables. 


Response variable Fixed effect Estimate Std. Error df t-value P-value 
SVL 0.378 0.296 59.929 1.276 0.207 
SVL 0.037 0.013 26.686 2.864 0.008 
SVL 0.044 0.012 61 3.724 <0.001 
SVL 0.019 0.007 60.859 2.550 0.013 
SVL 0.398 0.305 60.764 1.304 0.197 
SVL 0.042 0.012 61 3.371 0.001 
SVL 0.046 0.012 45.709 3.790 <0.001 
SVL 0.018 0.008 59.742 2:393 0.020 
SVL 0.5785 0.2946 60.272 1.963 0.054 
SVL 0.0408 0.0125 6l 3.277 0.002 
SVL 0.0477 0.0115 61 4.164 <0.001 
SVL 0.0222 0.0076 48.729 2.908 0.005 
SVL 0.32473 0.2908 59.73 1.117 0.269 
SVL 0.03857 0.0128 61 3.024 0.004 
SVL 0.04936 0.0117 61 4.214 <0.001 
SVL 0.01637 0.0071 60.178 2.316 0.024 
SVL 0.4464 0.274 60.989 1.629 0.109 
SVL 0.0348 0.012 61 2.824 0.006 
SVL 0.043 0.012 61 3.617 0.001 

0.0192 0.007 59.857 2.738 0.008 
0.44885 0.30192 60.948 1.487 0.142 
0.04008 0.0122 61 3.285 0.002 
0.04384 0.01159 61 3.782 <0.001 


0.02429 0.00718 16.36 3.385 0.004 


Continued Table S8 


Response variable Fixed effect Estimate Std. Error df t-value P-value 


SVL 0.3768 0.3037 60.56 1241 0.220 
SVL 0.03937 0.0123 61 3.201 0.002 
SVL 0.04196 0.01204 22.47 3.484 0.002 
SVL 0.02034 0.007622 57.48 2.668 0.010 


Note: V3d represents 3D volume of endocast; CL represents maximum endocast length; CW represents maximum endocast width; CH 
represents maximum endocast height; SVL represents body size. The bold font represents the significance (P < 0.05 or P < 0.001) of the 
correlation between endocast sizes and climate variables. 


